IRFs
Bernd Funovits
Summary
We obtain the IRFs for some of the best models (with respect to AIC/BIC/Shapiro-Wilk/Jarque-Bera/Ljung-Box). We use either estimated values for the static shock transmission matrix \(B\) or impose the long-run restriction suggested by Blanchard and Quah (1989).
In addition, we generate the IRFs of Blanchard and Quah (1989) and GMR.
Preliminaries
Load and plot data
tt = readRDS(paste0(params$PATH, params$JOBID, "/",
"e_residualcheck_end.rds")) %>% arrange(rk_bic)
DATASET = readRDS(params$PATH_DATASET)
dygraphs::dygraph(DATASET)DATASET = DATASET %>% as.matrix()
DIM_OUT = dim(DATASET)[2]
N_OBS = dim(DATASET)[1]Define some functions for IRFs
# Permute static shock transmission matrix B and change signs
permute_chgsign = function(irf_array,
perm = rep(1, dim(irf_array)[1]),
sign = rep(1, dim(irf_array)[1])){
dim_out = dim(irf_array)[1]
perm_mat = diag(dim_out)
perm_mat = perm_mat[,perm]
sign = diag(sign)
ll = map(1:dim(irf_array)[3], ~ irf_array[,,.x] %*% perm_mat %*% sign)
irf_array = array(0, dim = c(dim_out, dim_out, dim(irf_array)[3]))
for (ix in 1:dim(irf_array)[3]){
irf_array[,,ix] = ll[[ix]]
}
return(irf_array)
}
# Plot IRF
plot_irf = function(irf_array){
plot(pseries(irf_array %>% polm(), lag.max = 40))
}
# Calculate the cumulative sum of the IRF coefficients for a given element
irf_cumsum = function(irf_array, el_mat){
stopifnot("*el_mat* should be a matrix with two columns where each row contains the indices for which the cumulative sum should be calculated" =
dim(el_mat)[2] == 2)
for (ix in (1:nrow(el_mat))){
irf_array[el_mat[ix,1], el_mat[ix,2], ] = irf_array[el_mat[ix,1], el_mat[ix,2], ] %>% cumsum()
}
return(irf_array)
}
# Rotate the static shock transmission matrix such that there is no long-run effect
rotate_longrun = function(irf_array){
dim_out = dim(irf_array)[1]
stopifnot("Number of outputs must be equal to 2" = dim_out == 2)
k1 = freqresp(armamod(lmfd(polm(diag(dim_out)), polm(irf_array)), sigma_L = diag(dim_out)), n.f = 1)$frr
lq = lq_decomposition(matrix(c(Re(k1)), nrow = dim_out))
Q = lq$q
jj = map(1:dim(irf_array)[3], ~ irf_array[,,.x] %*% t(Q))
irf_array = array(0, dim = c(dim_out, dim_out,dim(irf_array)[3]))
for (ix in 1:dim(irf_array)[3]){
irf_array[,,ix] = jj[[ix]]
}
return(irf_array)
}Model statistics
tt %>%
select(-nr, -p_plus_q, -rk_mle, -contains("value"), -contains("cov"), -contains("pval")) %>%
filter(normality_flag == 0, lb_flag == 0) %>%
dtable()Blanchard Quah AR(8)
Function that calculates the cumulative IRF for some variables of choice. Here, we are interested in log(GNP) and not its change.
Estimate a VAR(8) model with the svars package.
vars_bq = VAR(DATASET, p = 8, type = "none")
vars_irf_bq = irf(vars_bq, n.ahead = 100)Transform it to rationalmatrices/RLDM style and create a first plot.
irf_bq = array(0, dim = c(2,2,101))
irf_bq[,1,] = t(vars_irf_bq$irf$rGDPgrowth_demeaned)
irf_bq[,2,] = t(vars_irf_bq$irf$unemp_detrended)
plot(pseries(polm(irf_bq), lag.max = 100))
Calculate the long-run response, LQ decompose it, rotate the static shock impact matrix, transform the transfer function s.t. the demand shock is transitory.
k1 = matrix(c(colSums(vars_irf_bq$irf$rGDPgrowth_demeaned),
colSums(vars_irf_bq$irf$unemp_detrended)),
nrow = 2)
lq = lq_decomposition(k1)
Q = lq$q
ii = map(1:101, ~ irf_bq[,,.x] %*% t(Q))
irf_bq_lr = array(0, dim = c(2,2,101))
for (ix in 1:101){
irf_bq_lr[,,ix] = ii[[ix]]
}
rm(k1, Q, ii, irf_bq)Check the result : transfer function evaluated at \(1\) should be triangular.
apply(irf_bq_lr, c(1,2), sum) %>% zapsmall()## [,1] [,2]
## [1,] 0.622691 0.000000
## [2,] -0.161077 5.558043Plot the rotated IRF
plot(pseries(irf_bq_lr %>% polm(), lag.max = 40))
irf_bq_lr_cumsum = irf_cumsum(irf_bq_lr, matrix(c(1,1,1,2), nrow = 2))Plot the rotated IRF where the response to (log) GNP is cumulated.
plot(pseries(irf_bq_lr_cumsum %>% polm(), lag.max = 40))
Permute columns and change sign such that the plots align with the ones in Blanchard and Quah.
jj = map(1:101, ~ irf_bq_lr_cumsum[,,.x] %*% matrix(c(0,-1,1,0), nrow = 2))
irf_bq_lr_cumsum_permuted = array(0, dim = c(2,2,101))
for (ix in 1:101){
irf_bq_lr_cumsum_permuted[,,ix] = jj[[ix]]
}
plot(pseries(irf_bq_lr_cumsum_permuted %>% polm(), lag.max = 40))
GMR VARMA(4,1, k = 1)
if(!file.exists("../local_data/gmr_results_bq.rds")){
load("~/r_projects/gmr_ssvarma/results/BQ/results.MLE.4lags.res")
write_rds(res.estim.MLE, "../local_data/gmr_results_bq.rds")
gmr_results_bq = res.estim.MLE
} else{
gmr_results_bq = read_rds("../local_data/gmr_results_bq.rds")
}Convert GMR results to rationalmatrices/RLDM style.
ar_gmr = gmr_results_bq$Phi
ma_gmr = gmr_results_bq$Theta
B_gmr = gmr_results_bq$C
polm_ar = array(c(c(diag(2)),
-c(ar_gmr)),
dim = c(2,2,5)) %>% polm()
polm_ma = array(c(c(diag(2)),
-c(ma_gmr)),
dim = c(2,2,2)) %>% polm()
B = gmr_results_bq$C
lmfd_obj = lmfd(polm_ar, polm_ma %r% B) Check poles and zeros
polm_ar %>% zeroes() %>% abs()## [1] 1.192445 1.192445 1.324240 1.447863 1.447863 2.216027 2.216027 4.102746polm_ma %>% zeroes() %>% abs()## [1] 0.565828 2.519183irf_gmr = pseries(lmfd_obj, lag.max = 40) %>% unclass()
plot_irf(irf_gmr)
irf_gmr_cumsum = irf_cumsum(irf_gmr,
el_mat = matrix(c(1,1,1,2), nrow = 2))
plot_irf(irf_gmr_cumsum)
Permute columns such that results align
irf_gmr_cumsum_permuted = permute_chgsign(irf_gmr_cumsum, perm = c(2,1))
plot_irf(irf_gmr_cumsum_permuted)
BF VARMA(1,2, kappa = 1, k = 0)
Extract appropriate inter-valued parameters.
tt %>%
filter(p==1 & q==2 & kappa==1 & k==0) %>%
select(p_plus_q, p,q, kappa, k, starts_with("rk_"), cov_el_sum, value_final, value_aic, value_bic, everything()) -> tt_bf_12_2 Create an ARMA-WHF model:
bf_params_12_2 = tt_bf_12_2$params_deep_final %>% .[[1]]
bf_tmpl_12_2 = tt_bf_12_2$tmpl %>% .[[1]]
write_rds(bf_params_12_2, "../local_data/p_whf/best_mod_rob/bf_params_12_2.rds")
write_rds(bf_tmpl_12_2, "../local_data/p_whf/best_mod_rob/bf_tmpl_12_2.rds")
bf_armawhfmod_12_2 = fill_tmpl_whf_rev(bf_params_12_2, bf_tmpl_12_2)Create IRF, cumulate GNP, and plot:
bf_irf_12_2 = irf_whf(bf_params_12_2, bf_tmpl_12_2, n_lags = 40) %>% unclass()
bf_irf_12_2_cumsum = bf_irf_12_2 %>% irf_cumsum(el_mat = matrix(c(1,1,1,2), nrow = 2))
bf_irf_12_2_cumsum %>% plot_irf()
bf_irf_12_2_cumsum_permuted = bf_irf_12_2_cumsumRotate such that BQ long-run restriction is satisfied, plot, sign-permute, and plot again.
bf_irf_12_2_lr = bf_irf_12_2 %>% rotate_longrun()
bf_irf_12_2_lr_cumsum = bf_irf_12_2_lr %>% irf_cumsum(el_mat = matrix(c(1,1,1,2), nrow = 2))
bf_irf_12_2_lr_cumsum %>% plot_irf()
bf_irf_12_2_lr_cumsum_permuted = bf_irf_12_2_lr_cumsum %>% permute_chgsign(perm = c(2,1),
sign = c(-1,1))
bf_irf_12_2_lr_cumsum_permuted %>% plot_irf()
BF VARMA(1,3, kappa = 1, k = 0)
Extract appropriate inter-valued parameters.
tt %>%
filter(p==1 & q==3 & kappa==1 & k==0) %>%
select(p_plus_q, p,q, kappa, k, starts_with("rk_"), cov_el_sum, value_final, value_aic, value_bic, everything()) -> tt_bf_13_2 Create an ARMA-WHF model:
bf_params_13_2 = tt_bf_13_2$params_deep_final %>% .[[1]]
bf_tmpl_13_2 = tt_bf_13_2$tmpl %>% .[[1]]
bf_armawhfmod_13_2 = fill_tmpl_whf_rev(bf_params_13_2, bf_tmpl_13_2)Create IRF, cumulate GNP, and plot:
bf_irf_13_2 = irf_whf(bf_params_13_2, bf_tmpl_13_2, n_lags = 40) %>% unclass()
bf_irf_13_2_cumsum = bf_irf_13_2 %>% irf_cumsum(el_mat = matrix(c(1,1,1,2), nrow = 2))
bf_irf_13_2_cumsum %>% plot_irf()
bf_irf_13_2_cumsum_permuted = bf_irf_13_2_cumsumRotate such that BQ long-run restriction is satisfied, plot, sign-permute, and plot again.
bf_irf_13_2_lr = bf_irf_13_2 %>% rotate_longrun()
bf_irf_13_2_lr_cumsum = bf_irf_13_2_lr %>% irf_cumsum(el_mat = matrix(c(1,1,1,2), nrow = 2))
bf_irf_13_2_lr_cumsum %>% plot_irf()
bf_irf_13_2_lr_cumsum_permuted = bf_irf_13_2_lr_cumsum %>% permute_chgsign(perm = c(2,1),
sign = c(-1,1))
bf_irf_13_2_lr_cumsum_permuted %>% plot_irf()
BF VARMA(3,1, kappa = 1, k = 0)
Extract appropriate inter-valued parameters.
tt %>%
filter(p==3 & q==1 & kappa==1 & k==0) %>%
select(p_plus_q, p,q, kappa, k, starts_with("rk_"), cov_el_sum, value_final, value_aic, value_bic, everything()) -> tt_bf_31_2 Create an ARMA-WHF model:
bf_params_31_2 = tt_bf_31_2$params_deep_final %>% .[[1]]
bf_tmpl_31_2 = tt_bf_31_2$tmpl %>% .[[1]]
bf_armawhfmod_31_2 = fill_tmpl_whf_rev(bf_params_31_2, bf_tmpl_31_2)Create IRF, cumulate GNP, and plot:
bf_irf_31_2 = irf_whf(bf_params_31_2, bf_tmpl_31_2, n_lags = 40) %>% unclass()
bf_irf_31_2_cumsum = bf_irf_31_2 %>% irf_cumsum(el_mat = matrix(c(1,1,1,2), nrow = 2))
bf_irf_31_2_cumsum %>% plot_irf()
Sign-permute such that shocks are labelled consistently.
bf_irf_31_2_cumsum_permuted = bf_irf_31_2_cumsum %>% permute_chgsign(perm = c(2,1),
sign = c(1,1))
bf_irf_31_2_cumsum_permuted %>% plot_irf()
Rotate such that BQ long-run restriction is satisfied, plot, sign-permute, and plot again.
bf_irf_31_2_lr = bf_irf_31_2 %>% rotate_longrun()
bf_irf_31_2_lr_cumsum = bf_irf_31_2_lr %>% irf_cumsum(el_mat = matrix(c(1,1,1,2), nrow = 2))
bf_irf_31_2_lr_cumsum %>% plot_irf()
bf_irf_31_2_lr_cumsum_permuted = bf_irf_31_2_lr_cumsum %>% permute_chgsign(perm = c(2,1),
sign = c(-1,1))
bf_irf_31_2_lr_cumsum_permuted %>% plot_irf()
BF VARMA(1,1, kappa = 1)
Extract appropriate inter-valued parameters.
tt %>%
filter(p==1 & q==1 & kappa==1 & k==0) %>%
select(p_plus_q, p,q, kappa, k, starts_with("rk_"), cov_el_sum, value_final, value_aic, value_bic, everything()) -> tt_bf_11_2 Create an ARMA-WHF model:
bf_params_11_2 = tt_bf_11_2$params_deep_final %>% .[[1]]
bf_tmpl_11_2 = tt_bf_11_2$tmpl %>% .[[1]]
bf_armawhfmod_11_2 = fill_tmpl_whf_rev(bf_params_11_2, bf_tmpl_11_2)Create IRF, cumulate GNP, and plot:
bf_irf_11_2 = irf_whf(bf_params_11_2, bf_tmpl_11_2, n_lags = 40) %>% unclass()
bf_irf_11_2_cumsum = bf_irf_11_2 %>% irf_cumsum(el_mat = matrix(c(1,1,1,2), nrow = 2))
bf_irf_11_2_cumsum %>% plot_irf()
Sign-permute such that shocks are labelled consistently.
bf_irf_11_2_cumsum_permuted = bf_irf_11_2_cumsum %>% permute_chgsign(perm = c(2,1),
sign = c(-1,-1))
bf_irf_11_2_cumsum_permuted %>% plot_irf()
Rotate such that BQ long-run restriction is satisfied, plot, sign-permute, and plot again.
bf_irf_11_2_lr = bf_irf_11_2 %>% rotate_longrun()
bf_irf_11_2_lr_cumsum = bf_irf_11_2_lr %>% irf_cumsum(el_mat = matrix(c(1,1,1,2), nrow = 2))
bf_irf_11_2_lr_cumsum %>% plot_irf()
bf_irf_11_2_lr_cumsum_permuted = bf_irf_11_2_lr_cumsum %>% permute_chgsign(perm = c(2,1),
sign = c(-1,1))
bf_irf_11_2_lr_cumsum_permuted %>% plot_irf()
BF VARMA(4,1, kappa = 0, k = 1): GMR’s model
Extract appropriate inter-valued parameters.
tt %>%
filter(p==4 & q==1 & kappa==0 & k==1) %>%
select(p_plus_q, p,q, kappa, k, starts_with("rk_"), cov_el_sum, value_final, value_aic, value_bic, everything()) -> tt_bf_41_1 Create an ARMA-WHF model:
bf_params_41_1 = tt_bf_41_1$params_deep_final %>% .[[1]]
bf_tmpl_41_1 = tt_bf_41_1$tmpl %>% .[[1]]
bf_armawhfmod_41_1 = fill_tmpl_whf_rev(bf_params_41_1, bf_tmpl_41_1)Create IRF, cumulate GNP, and plot:
bf_irf_41_1 = irf_whf(bf_params_41_1, bf_tmpl_41_1, n_lags = 40) %>% unclass()
bf_irf_41_1_cumsum = bf_irf_41_1 %>% irf_cumsum(el_mat = matrix(c(1,1,1,2), nrow = 2))
bf_irf_41_1_cumsum %>% plot_irf()
Sign-permute such that shocks are labelled consistently.
bf_irf_41_1_cumsum_permuted = bf_irf_41_1_cumsum %>% permute_chgsign(perm = c(2,1))
bf_irf_41_1_cumsum_permuted %>% plot_irf()
Rotate such that BQ long-run restriction is satisfied, plot, sign-permute, and plot again.
bf_irf_41_1_lr = bf_irf_41_1 %>% rotate_longrun()
bf_irf_41_1_lr_cumsum = bf_irf_41_1_lr %>% irf_cumsum(el_mat = matrix(c(1,1,1,2), nrow = 2))
bf_irf_41_1_lr_cumsum %>% plot_irf()
bf_irf_41_1_lr_cumsum_permuted = bf_irf_41_1_lr_cumsum %>% permute_chgsign(perm = c(2,1),
sign = c(-1,1))
bf_irf_41_1_lr_cumsum_permuted %>% plot_irf()
Summary
list_mods = tibble(BQ_AR8 = irf_bq_lr_cumsum_permuted[,,1:41] %>% list(),
GMR_ARMA41_1 = irf_gmr_cumsum_permuted %>% list(),
BF_ARMA11_2 = bf_irf_11_2_cumsum_permuted %>% list(),
BF_ARMA11_2_lr = bf_irf_11_2_lr_cumsum_permuted %>% list(),
BF_ARMA12_2 = bf_irf_12_2_cumsum_permuted %>% list(),
BF_ARMA12_2_lr = bf_irf_12_2_lr_cumsum_permuted %>% list(),
BF_ARMA13_2 = bf_irf_13_2_cumsum_permuted %>% list(),
BF_ARMA13_2_lr = bf_irf_13_2_lr_cumsum_permuted %>% list(),
BF_ARMA31_2 = bf_irf_31_2_cumsum_permuted %>% list(),
BF_ARMA31_2_lr = bf_irf_31_2_lr_cumsum_permuted %>% list(),
BF_ARMA41_1 = bf_irf_41_1_cumsum_permuted %>% list(),
BF_ARMA41_1_lr = bf_irf_41_1_lr_cumsum_permuted %>% list())
list_mods %>% pivot_longer(everything()) %>%
mutate(Demand2GNP = map(value, ~.x[1,1,]),
Demand2Unempl = map(value, ~.x[2,1,]),
Supply2GNP = map(value, ~.x[1,2,]),
Supply2Unempl = map(value, ~.x[2,2,])) %>%
select(-value) %>%
pivot_longer(c("Demand2GNP", "Demand2Unempl", "Supply2GNP", "Supply2Unempl"),
names_to = "Response_Type", values_to = "Impact") %>%
separate(Response_Type, into = c("Shock", "Output"), sep = "2") %>%
unnest_longer(Impact, indices_include = TRUE) %>%
rename(Lag = Impact_id,
Model = name) %>%
mutate(Lag = Lag-1) -> tibble_irftibble_irf %>%
# filter(Model != "BF_ARMA11_2") %>%
ggplot(aes(x = Lag, y = Impact, color = Model)) +
geom_line() + geom_point() + facet_grid(Output~Shock) -> ply
plotly::ggplotly(ply)knitr::knit_exit()